function [p_value] = SongLeungU(X,Y,R,b,ind,h0)
% [rej,stat] = SongLeungU(X,Y,R,b,ind,h0,alpha) returns Song/Leung
% test results from the "U-statistic" version of their procedure
% for the coefficients corresponding to the entries of ind and null
% hypothesis specified by h0.
% Looking only at marginal inference on each element in ind, not joint
% hypothesis testing.
% X - design matrix
% Y - outcome
% R - number of subsamples
% b - subsample size
% ind - coefficients for which a confidence interval is desired
% h0 - hypothesized value for parameters of interest
% alpha - level of test
% seed - for replicability, might want to set seed of random number
% generator

if nargin > 7 
    if ~isempty(seed)
        rng(seed);
    end
end

[n,p] = size(X);

pseudoobs = n*(X.*(Y*ones(1,p)))/(X'*X);  

S = ((pseudoobs - ones(n,1)*mean(pseudoobs))'*(pseudoobs - ones(n,1)*mean(pseudoobs)))/n;

p0 = numel(ind);
stat = zeros(p0,1);
for jj = 1:ceil(R)
    samplejj = randsample(n,ceil(b));
    for kk = 1:p0
        pseudoobsjjkk = pseudoobs(samplejj,kk)-h0(kk);
        stat(kk) = stat(kk) + (sum(pseudoobsjjkk)^2 - sum(pseudoobsjjkk.^2))/S(kk,kk);
    end
end
stat = stat/(sqrt(2*ceil(R))*ceil(b));

p_value = 2*(1-normcdf(abs(stat)));
% rej = stat > norminv(1-alpha);
            
    
        